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ABSTRACT 


Few  previous  investigators  have  attempted  to  experimentally 
study  unsteady  state  flow  of  gas  through  porous  media.  Instead, 
they  have  predicted  transient  behaviour  by  various  approximations 
to  the  non-linear  partial  differential  equation  describing  such  flow. 

In  the  course  of  this  investigation,  the  transient  behaviour 
was  measured  directly  by  flowing  nitrogen  through  relatively  long 
full-sized  core  samples.  The  experimental  results  were  scruti¬ 
nized  and  compared  with  numerical  solutions  to  a  mathematical 
model.  The  boundary  conditions  simulated  were: 

1.  Constant  terminal  pressure  with  constant  pressure 

at  the  external  boundaiy  . 

2 .  Constant  terminal  pressure  with  a  sealed  external 
boundary. 

The  results  obtained  indicate  the  shortcomings  of  certain 
types  of  flow  systems,  the  desirability  of  long  flow  systems  and 
the  behaviour  of  pressure  transients  in  both  homogeneous  and  hetero¬ 
geneous  systems.  Furthermore,  the  investigation  illustrates  the 
approximation  to  experimental  results  obtained  by  means  of  the 


mathematical  model. 
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INTRODUCTION 

Petroleum  industry  studies  involving  unsteady  state  or 
transient  flow  generally  are  concerned  with  analysis  of  well  tests 
or  water  influx.  The  well  tests  include  drawdown,  buildup,  inter¬ 
ference  and  back-pressure  tests  which  can  be  utilized  to  evaluate 
such  quantitative  information  about  a  well  and  reservoir  as  well¬ 
bore  damage,  effective  permeability  within  the  drainage  radius, 
static  reservoir  pressure,  distances  to  boundaries,  inter-well 
communication,  reservoir  size  and  capabilities.  These  evaluations 
are  simply  applications  of  solutions  of  modified  forms  of  the  con¬ 
tinuity  equation.  A  very  general  form  of  the  continuity  equation  is 
given  below. 

-  -t  _3£ 

at 

In  the  case  of  water  influx,  the  aquifer  response  is  related 

to  the  variations  of  the  average  boundary  pressure  by  solutions  to 

(31) 

the  continuity  equation.  Van  Everdingen  and  Hurst  have  solved 
this  form  of  the  continuity  equation  for  influx  and  have  presented 
the  results  in  convenient  tabular  and  graphical  forms  in  terms  of 
dimensionless  quantities. 

The  objective  of  this  study  was  to  experimentally  measure 
and  examine  unsteady  state  or  transient  gas  flow  in  porous  media. 


pv 


I 
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In  general  terms  this  type  of  fluid  flow  can  be  described  by  the 
continuity  equation.  Substitution  of  an  equation  of  motion  and 
an  equation  of  state  into  the  continuity  equation  yields  a  form  of 
the  diffisivity  equation  utilized  in  this  study.  This  form  of  the 
diffusivity  equation  describes  isothermal,  viscous  flow  of  a  real 
gas  in  a  linear,  horizontal  system, 

Darcy’s  law  which  was  used  as  the  equation  of  motion  is 
limited  to  the  viscous  flow  regime  and  does  not  account  for  slip 
and  inertial  energy  losses.  The  seriousness  of  this  simplification 
was  given  a  cursory  examination  by  comparing  the  numerical 
results  obtained  from  a  mathematical  model  and  the  experimental 
results.  The  model  incorporates  the  variation  of  viscosity  and 
compressibility  with  pressure  and  the  variation  of  absolute  per¬ 


meability  with  position. 
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THEORY  AND  LITERATURE  REVIEW 

The  movement  of  gas  through  porous  media  has  been 
extensively  investigated,  and  from  these  investigations  three  flow 
regions  have  been  established  for  separate  consideration:  (a) 
molecular  streaming,  (b)  viscous,  and  (c)  turbulent  flow. 

Scheidegger  summarizes  significant  theoretical  and 
experimental  contributions  in  these  vital  regions,  These  early 
examinations  deal  with  steady  state  flow  primarily,  wherein  the 
pressure  and  fluid  velocity  at  every  point  throughout  the  system 
adjust  instantaneously  to  either  a  pressure  or  flow  rate  change 
in  some  other  part  of  the  system. 

More  recently,  the  unsteady  state  or  transient  flow  of  gas 
has  been  receiving  considerable  attention,  especially  in  relation 
to  back  pressure  testing  and  pressure  build  up  analysis.  Expansion 
of  compressed  reservoir  fluids  during  production  results  in  un¬ 
steady  state  flow.  When  a  well  is  opened  for  flow,  a  pressure 
drop  occurs  at  the  wellbore  causing  fluid  in  the  vicinity  of  the  well¬ 
bore  to  expand  and  flow  towards  it.  As  a  result  of  this  fluid  loss, 
the  reservoir  pressure  more  remote  from  the  wellbore  declines, 
permitting  the  fluid  there  to  expand  and  to  flow  towards  the  well. 
Consequently  the  drainage  area  continues  to  expand,  and  the 
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velocity  and  pressure  are  continuously  changing  at  all  points. 
Hence,  in  an  elemental  volume  of  porous  media,  the  flow  rate 
into  the  media  is  not  the  same  as  that  out  of  the  media. 

An  equation  which  properly  describes  transient  flow  must 
firstly  involve  a  material  balance: 


Amount  of 

— 

Amount  of 

— 

Increase  of  Mass 

Mass  Input 

Mass  Output_ 

Within  the  Region 

Based  on  the  conservation  of  mass  for  isothermal  fluid 
flow  through  a  porous  medium,  the  material  balance  takes  the 
form  of  the  well-known  continuity  equation: 

=  -  4>  J P_  (i) 

dt 

The  vector  in  Equation  1,  which  represents  the  volu¬ 
metric  rate  of  flow  per  unit  cross-sectional  area,  can  be 
expressed  by  Darcy’s  Law  if  the  flow  is  viscous: 

v  =  -  k(p)  P  V  l  (2) 

»(P) 

n  7) 

Hubbert  has  studied  Darcy’s  Law  and  its  implications 
and  has  indicated  that: 

$  =  i  &  +  gz 

o  r 


(3) 
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Generally  gravitational  effects  are  ignored,  and  Darcy’s 
Law  becomes : 


v  =  _  k(p)  V  p  (4) 

u(p) 

At  high  velocities,  Darcy’s  Law  no  longer  expresses  the 
relationship  between  the  flow  rate  and  the  pressure  gradient,  but 
must  be  replaced  by  an  equation  incorporating  a  quadratic  velocity 
term: 

Dv  Ivl  +  t  =  -  k(P)  VP  (5> 

u(p) 


However,  assuming  viscous  flow,  the  combination  of  the 
continuity  equation  (Equation  1)  and  Darcy’s  Law  (Equation  4) 
yields  a  differential  equation  describing  fluid  flow  in  a  porous 
medium. 


V 


gjsfri 

_  u(p) 


VP 


(6) 


The  final  form  of  the  equation  which  adequately  describes 
unsteady  state  flow  through  porous  media  depends  on  the  fluid  and 
the  equation  of  state.  For  real  gases: 


<°g 


_M _ p_ 

RT  z  (p) 


(7) 


■ 
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Eliminating  density  (p)  from  Equation  6  and  assuming 
isothermal  flow  will  yield; 

V*  fk(P) 

Lp0(p) 25 (p) 

Equation  8  with  necessary  boundary  and  initial  conditions 
along  with  appropriate  correlations  relating  permeability,  viscosity, 
and  compressibility  factor  to  pressure,  will  mathematically  describe 
the  flow  of  gas  through  a  porous  medium.  The  proposed  experimental 
study  involves  no  more  than  analysing  the  pressure  behaviour  as 
a  function  of  both  position  and  time  and  comparing  this  with  numerical 
solutions  to  Eq,  8, 

A  search  through  the  literature  revealed  that  the  experimen¬ 
tal  examination  of  unsteady  state  gas  flow  is  very  limited.  Therefore, 
despite  this  being  primarily  an  experimental  study,  a  review  of  the 
theoretical  investigations  related  to  this  type  of  flow  was  undertaken 
to  better  understand  the  implications  of  Equation  (8)  and  its  associated 
assumptions.  Most  of  the  theoretical  solutions  involved  the  simplifying 
assumptions  of  constant  (a)  permeability,  (b)  viscosity  and  (c)  compres¬ 
sibility  factor,  which  reduces  Equation  (8)  to: 


PVP 


=  P 


-pq 

t  Lz(p)j 


(8) 


P?P 


ji  p 
k  d  t 


(9) 
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The  methods  included  pseudo-steady  state,  analytical, 

graphical,  numerical  and  analogue  solutions,  A  review  of 

graphical  and  analogue  solutions  will  not  be  presented  here,  but 

can  be  found  in  Katz  et  al  One  of  the  earliest  attempts  to 

solve  Equation  (9)  was  made  by  Heatherington,  MacRoberts  and 
(15) 

Huntington  .  Realizing  that  nearly  all  phenomena  of  practical 

interest  in  gas  production  involved  unsteady  state  flow,  and  that 

the  pertinent  equation  was  not  mathematically  reducible  to  simple 

expressions  void  of  infinite  series,  Heatherington  et  al  proposed  an 

approximate  solution  based  upon  "infinite  series  of  steady  states" . 

With  time  no  longer  an  independent  variable,  the  equation  was 

solved  by  simple,  explicit  integration.  The  authors  were  able 

to  conclude  from  a  comparison  of  experimental  and  theoretical 

data  that  the  general  back -pressure  equation  was  reasonably 

accurate,  provided  the  initial  transient  period  was  not  too  long. 

(24) 

Muskat  similarly  obtained  solutions  to  Equation  (9) 
by  assuming  a  "succession  of  steady  states".  This  implies  that 
the  velocity  of  the  propogation  of  disturbances  in  a  porous  medium 
is  infinite,  while  the  mass  fluid  content  variation  of  the  system  is 
negligible  compared  to  the  steady  state  flux  through  the  system. 


-.<fU  ( 1  i  juifi  c*  s  f  <»J  hyjfcqmoo  al  Uy'tgM 
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However,  as  indicated  by  Muskat,  these  two  assumptions  could  be 
most  inaccurate  if  the  system  were  very  large  or  where  abnormally 
high  compressibilities  were  encountered.  He  proceeds  to 
illustrate  by  example  how  an  oil  system  containing  as  little  as  4,5 
per  cent  free  gas  would  have  sufficiently  high  compressibility  to 
make  the  application  of  "succession  of  steady  states”  invalid. 

Utilizing  the  steady  state  formula  for  radial  gas  flow, 

(3) 

Aronofsky  and  Jenkins  '  approximated  the  transient  solution  to 
Equation  (9)  by  suitably  defining  an  effective  drainage  radius. 
Furthermore,  these  authors  were  able  to  conclude  that  the  solu¬ 
tions  for  the  liquid  flow  case  could  be  used  to  generate  approximate 

solutions  for  the  constant  rate  of  production  of  ideal  gases.  Extending 

(1) 

this  approach,  Al-Hassainy,  Ramey  and  Crawford  approximated 
Equation  (8)  by  introducing  a  variable  change  applicable  to  the 
flow  of  real  gases  through  porous  media.  This  change  was  used 
as  a  pseudo-pressure  consisting  of  pressure  as  well  as  the  pressure 
dependent  quantities  of  viscosity  and  compressibility  factor.  The 
final  form  of  the  diffusivity  equation  simply  has  pressure  squared 
replaced  by  the  pseudo-pressure  term. 


■ 
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A  summary  of  the  more  significant  analytical  solutions  is 

(18) 

presented  by  Katz  et  al  in  tabular  and  graphical  form  for  both 

linear  and  radial  flow.  The  mathematical  derivation  for  linear  flow 

(9) 

was  obtained  from  heat  flow  studies  by  Churchill  ,  while  the 

radial  flow  solutions  were  developed  by  Van  Everdingen  and  Hurst^^, 
(28) 

Roberts  obtained  an  approximate  analytical  solution  to  the  one 
dimensional  form  of  Equation  (9),  He  employed  a  stepwise  forward 
integration  with  respect  to  time  to  solve  a  linearized  form  of  the 
partial  differential  equation.  By  comparing  the  results  to  those 
obtained  when  assuming  the  gas  obeyed  the  linear  heat  conduction 
equation,  he  showed  that  the  non-linearity  reduced  the  rate  of 
pressure  decline,  Kidder  applying  perturbation  techniques 
and  utilizing  the  well-known  Boltzmann  transformation,  presented 
an  analytical  solution  for  gas  flow  in  a  semi-infirtite  porous  medium. 
With  the  advent  of  digital  computers,  numerical  solutions 
using  both  finite  difference  equations  and  computing  techniques 
have  become  very  popular,  Aronofsky  and  Jenkins  ^  obtained 
solutions  to  the  one  dimensional  form  of  Equation  (9)  by  employing 
a  forward  difference  explicit  (FDE)  approximation  modified  with 
respect  to  time  to  insure  stability.  They  observed  from  finite 


. 
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tube  studies  that  the  non-linear  character  of  the  equation  caused  the 
solution  to  approach  steady  state  more  rapidly  for  the  larger  ratios 
of  initial  pressure  to  constant  terminal  pressure.  Similarly,  Bruce, 
Peaceman,  Rachford  and  Rice  utilizing  FDE,  obtained  solutions 
for  linear  flow.  However,  they  also  determined  by  an  implicit 
technique,  the  solutions  for  linear  and  radial  flow. 

The  linear  and  radial  forms  of  Equation  (8),  with  perme¬ 
ability  constant,  were  solved  by  Aronsfsky  and  Ferris  ^  and 

(5) 

by  Aronsfsky  and  Porter  respectively.  The  radial  consideration 
was  modified  to  insure  stability  by  assuming  a  steady  state  core 
around  the  wellbore.  These  authors  observed  that  the  influence 
of  varying  gas  properties  had  a  significant  effect  on  the  time 
dependency  of  the  pressure. 

(25,  26) 

Quon,  Dranchuk,  Allada  and  Leung  discussed 

five  numerical  techniques  for  handling  Equation  (8)  in  a  two 
dimensional  system.  These  techniques  are:  (a)  forward  difference 
explicit  (FDE),  (b)  backward  difference  implicit  (BDI),  (c)  Crank- 
Nicholson  implicit  (CNI),  (d)  alternating  -  direction  explicit  pro¬ 
cedure  (ADEP)  and  (e)  alternating-direction  implicit  procedure 
(ADIP).  Quon  et  al  concluded  ADEP  offers  computer  efficiency 


over  the  implicit  procedures  and  has  a  stability  advantage  over  FDE. 
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(g) 

Carter  examined  three  methods  of  solving  Equation  (8) 
in  a  two  dimensional  model  and  observed  that  implicit  procedures 
of  the  Saul’ev  type  are  stable  for  practical  time  step  sizes.  Further¬ 
more  from  a  computing  time  viewpoint,  these  procedures  are  more 

desirable  than  the  explicit  techniques. 

(13) 

Eilerts  presents  solutions  to  Equation  (8)  for  flow  of 
gas  condensate  fluids  in  linear  and  radial  systems.  Viscosity 
and  compressibility  are  treated  as  functions  of  pressure  while 
permeability  is  a  function  of  pressure  and  position.  The  resulting 
partial  differential  equation  is  written  in  the  form  of  finite  differences 
and  solved  by  available  explicit  and  implicit  procedures. 

This  study  involves  a  comparison  of  experimental  unsteady- 
state  flow  data  with  that  predicted  by  a  one -dimensional  mathematical 
model.  Since  it  was  desirable  that  the  model  incorporate  variability 
of  viscosity  and  compressibility  with  pressure  and  the  variability 
of  permeability  with  distance,  the  diffusivity  equation  (Eq,  8)  was 
modified  to  give: 


3 

r  k(x)  P  a  pi 

-  $  d 

P 

3x 

_^(P)  z(p)  3x__ 

9t 

_Z(P)_ 
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Variation  of  viscosity  and  compressibility  with  pressure  was  provided 
by  correlations  found  in  the  literature  while  position  variation 

of  permeability  was  obtained  from  steady  state  calculations. 

o 

Furthermore  the  equation  was  linearized  with  respect  to  p  by  intro¬ 
ducing  a  mean  pressure: 


01) 

Finally,  specification  of  the  initial  and  boundary  conditions  to  be 
imposed  for  the  constant  terminal  pressure  case  completes  the 
model. 

Constant  Pressure  at  External  Boundary 
Initial  Condition 

p(x,  0)  =  p.  os?x=sL 

Boundary  Conditions 

P(o,t)  =pw  b>o 

p(L,t)  =  p.  t>o 

Sealed  External  Boundary 

Initial  Condition 

p(x,  o)  =  p. 


a 

K(x)  Qv* 

=  4 

a 

2~ 

P 

8x 

^(P)  z(P)  3x 

P 

3t 

_z(P) 

o^x^L 
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Boundary  Conditions 

P(o,  t)  =  pw 

t>o 

=  0 

t»o 

dx  J 

x=L 

Crank-Nicholson  finite  difference  equations  as  proposed  by 
(13) 

Eilerts  v  and  the  square  root  algorithm  extended  to  non-symmetric 
matrices,  were  utilized  in  the  solution  to  the  diffusivity  equation. 

Appendix  A  presents  the  derivation  of  the  PDE  while 
Appendix  B  illustrates  the  details  of  the  discretization  of  the 
finite  difference  equations  which  were  solved  on  an  IBM  360/67 
digital  computer.  Details  of  the  program  are  described  and 


presented  in  Appendix  C. 


. 
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EXPERIMENTAL  APPARATUS 


In  conjunction  with  the  measurement  of  pressure  transients, 
the  experimental  determination  of  the  absolute  permeability  and 
porosity  necessitated  measurements  of  pressure,  temperature  and 
flow  rate.  To  accomplish  this  a  wide  variety  of  equipment  was 
required. 

A  porosimeter  was  constructed  to  facilitate  pore  volume 
measurement  of  the  large,  encased  core  samples.  The  essential 
components  of  this  equipment  as  shown  in  Figure  I  are  a  calibrated 
Marsh  gauge  to  measure  the  initial  pressure  of  the  system,  an 
expansion  chamber,  a  mercury  manometer  to  measure  the  expanded 
pressure  and  a  vacuum  pump. 

Figure  II  is  a  composite  diagram  of  the  experimental  equip¬ 
ment  used  in  procuring  permeability  and  transient  data.  This  is 

(7) 

essentially  the  same  experimental  apparatus  employed  by  Cachero 
with  a  few  modifications.  These  changes  include  (1)  a  Moore  Nullmatic 
pressure  regulator  (2)  a  sixth  transducer  (3)  one  half  inch  O.  D.  flow 
line  on  the  downstream  side  of  the  core.  Constant  boundary  pressures 
were  ensured  by  the  first  and  third  modifications  while  the  second 
provided  a  check  of  pressure  drop  across  the  core  as  measured 


by  the  pressure  gauges. 


. 

« 
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DIAGRAM  OF  EXPERIMENTAL  APPARATUS 
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A  cylinder  charged  with  commercially  pure  nitrogen  was 
equipped  with  a  regulator  to  provide  gas  at  a  relatively  constant 
pressure.  Additional  control  was  supplied  by  a  pressure  cylinder 
functioning  as  a  surge  tank  and  the  nullmatic  pressure  regulator 
installed  upstream  from  the  core  holder.  Preceding  this  regulator 
and  immediately  downstream  from  the  core  holder  were  flowline 
iron-constantan  thermocouples  for  temperature  measurement. 
Bourdon  tube  pressure  gauges  tapped  into  each  end  plate  of  the 
core  holder  provided  means  of  measuring  pressure.  Flow  rate 
measurements  were  made  by  a  three-range  Voi-O-Flo  meter  and 
a  Precision  Instrument  wet  test  meter. 

Finally,  the  pressure  distribution  along  the  length  of  the 
core  was  obtained  using  calibrated  high  pressure  Statham  differ¬ 
ential  pressure  transducers  whose  output  voltage  was  continuously 
recorded  by  a  multi-channel  Honeywell  906  A  visicorder,  Input 
voltage  to  the  transducers  was  provided  by  a  constant  voltage 
source  and  checked  with  a  digital  voltmeter.  With  the  exception 
of  the  nitrogen  cylinder,  pressure  gauges  and  the  flow  meter,  the 
apparatus  was  enclosed  within  a  thermostatically  controlled  cabinet. 
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EXPERIMENTAL  PROCEDURE 


CORE  MOUNTING  AND  CLEANING  PROCEDURE 

A  berea  sandstone  model  labelled  Berea  1  in  this  study  and 

a  Series  core  were  those  prepared  by  Cachero.  In  addition  to  the 

aforementioned  cores,  full  size  Berea  sandstone  and  limestone 

samples  were  selected  for  non-steady  state  study.  These  cores 

(14) 

were  sheared  off  at  a  desired  length  as  described  by  Hamilton 
in  an  attempt  to  avoid  possible  mechnical  end  effects.  After  being 
centered  in  casing  by  a  mounting  board  and  wedges,  the  cores  were 

(23) 

set  using  the  identical  technique  and  epoxy  resin  employed  by  Mackett 
These  two  cores  plus  the  Berea  1  core  were  set  in  five  and  one-half  inch 
casing.  Thus  the  one  set  of  end  plates  would  serve  all  three  cores. 

The  Series  core  was  mounted  in  four  and  one-half  inch  casing  and 
was  equipped  with  a  set  of  screw-on  end  plates.  All  cores  were 
tapped  with  five  equally  spaced  taps  drilled  through  the  casing  and 
epoxy  to  the  core  surface.  The  Berea  1  and  Series  core  are 
illustrated  in  Figure  3. 

The  mounted  cores  were  flushed  with  isoctane,  chased  by 
propane  and  dried  with  nitrogen  prior  to  testing. 


*  i 


■ 


SERIES  CORE  AND  HOLDER 
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EFFECTIVE  POROSITY  DETERMINATION 

Bulk  volume  determination  of  the  mounted  cores  was  per¬ 
formed  by  direct  measurement  of  diameter  and  length,  while  the 
gas  expansion  porosimeter  provided  the  gross  grain  volume.  Four 
trials,  each  at  a  different  initial  pressure  ensured  that  the  pore 
volume  determination  was  not  pressure  dependent.  That  is, 
expansion  of  the  casing  and/or  epoxy  annulus  while  the  system 
was  pressurized  did  not  significantly  affect  porosity. 

FLOW  TESTS 

Flow  data  were  obtained  on  all  cores  using  the  experimental 
equipment  illustrated  in  Figure  2„  Prior  to  testing,  the  system 
was  checked  under  pressure  for  leaks  using  a  soap  solution  and 
observing  a  possible  decline  in  pressure  gauge  readings.  During 
permeability  flow  tests  the  tranducers  were  isolated  from  the 
system  to  avoid  excessive  pressure  differentials  on  any  one 
transducer, 

A  permeability  determining  flow  test  commenced  with  the 
setting  of  the  nullmatic  pressure  regulator  at  the  desired  upstream 
pressure.  Periodic  observation  of  the  upstream  and  downstream 
pressure  and  the  flow  rate  indicated  when  steady  state  flow  conditions 
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had  been  attained  and  various  observations  could  be  recorded.  At 
pressures  in  excess  of  the  operating  limits  of  the  nullmatic  regulator, 
it  was  removed  with  upstream  pressures  being  subsequently  controlled 
by  the  regulator  on  the  nitrogen  bottle. 

Each  experimental  run  included  the  following  recordings 
(1)  upstream  core-face  pressure,  (2)  downstream  core-face  pressure, 
(3)  upstream  thermocouple  reading,  (4)  downstream  thermocouple 
reading,  (5)  flow  rate  reading,  (6)  barometric  pressure,  and  (7) 
room  temperature  reading. 

The  unsteady  state  flow  tests  were  conducted  using  the 
same  experimental  setup;  however,  in  this  case  the  transducers 
were  in  communication  with  the  pressure  taps  but  isolated  from 
each  other. 

It  was  desirable  to  examine  two  different  boundary  conditions. 
Firstly,  constant  terminal  pressure  with  a  constant  pressure  at 
the  external  boundary  and  secondly,  constant  terminal  pressure  with 
a  sealed  external  boundary.  The  first  boundary  condition  was 
simulated  by  closing  the  downstream  toggle  valve  and  adjusting  the 
nullmatic  pressure  regulator  to  pressurize  and  provide  the  system  with 


the  desired  initial  pressure.  Then,  the  visicorderfs  chart  time  and 


. 
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speed  were  set  and  the  downstream. snap  action  valve  was  opened. 
Closure  of  the  valve  immediately  upstream  from  the  core  prior  to 
drawdown  provided  data  for  the  second  set  of  boundary  conditions. 

Typical  experimental  data  recorded  during  each  test 
immediately  prior  to  the  testveae  (1)  visicorder  time  and  speed 
setting,  (2)  upstream  core-face  pressure,  (3)  downstream  core¬ 
face  pressure,  (4)  upstream  thermocouple  reading.  (5)  downstream 
thermocouple  reading,  (6)  barometric  pressure,  and  (7)  room 
temperature  reading.  Furthermore  the  upstream  core-face 
pressure  was  recorded  at  various  times  to  provide  a  chedfer  and 
the  upstream  and  downstream  core-face  pressures  were  again 
recorded  once  steady  state  conditions  had  been  achieved. 

Initially,  flow  rates  at  various  times  during  each  test 
were  to  be  recorded;  however,  it  was  found  that  with  the  additional 
equipment  downstream  of  the  core-face,  constant  terminal  pressure 
could  not  be  achieved  prior  to  the  transients  reaching  the  upstream 


core  face. 


' 
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PROCESSING  OF  EXPERIMENTAL  RESULTS 

POROSITY  MEASUREMENTS 

Table  I  was  prepared  from  direct  measurements  and  from 
evaluated  porosimeter  data. 

TABLE  1 


Core 

Diameter 

(Inches) 

Length 

(Inches) 

Bulk  Volume  Porosity 
(Cu.  In.)  (Fraction) 

Berea  1 

3.51 

34.3 

320.8 

.208 

Berea  1A 

3.51 

20.8 

201.3 

.196 

Berea  2 

3.51 

21.4 

209.7 

.237 

Limestone 

3.51 

21.5 

203.1 

.213 

Series 

1.997 

54.0 

169.1 

.264 

Cores  Berea  1  and  Berea  1A  are  the  same  test  samples.  The 
Berea  1  sample  was  that  studied  by  Caehero  but  with  three-eights  of 
one  inch  fractured  off  each  end  to  examine  mechanical  end  effects. 
Berea  1A  was  the  same  sample  with  an  additional  six  and  one -half 
inches  removed  from  each  end. 

STEADY  STATE  MEASUREMENTS 

The  liquid  or  absolute  permeability  of  each  core  sample 
was  initially  determined  by  applying  the  appropriate  form  of  Darcy's 


.  ,  :  j  , ' .  jj  v.-i  L  .'Hi  -  teb  x  xWl  mw 
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law  to  test  data  in  the  viscous  flow  region  and  correcting  for 
slippage  as  suggested  by  Klinkenberg^21^  A  computer  program 
eliminated  most  of  the  manual  calculations. 

Upon  examining  the  back -pressure  plots  corrected  for 

(12) 

slip  as  per  Dranchuk  and  Kolada  ,  it  was  evident  that  some 
viscous  inertial  data  was  being  included  in  the  original  permeability 
determinations  for  the  Limestone  sample.  Furthermore,  flow  test 
data  for  the  Series  core  were  almost  entirely  in  the  visco-inertial 
region.  The  former  problem  was  corrected  by  eliminating  the 
visco-inertial  data  from  the  calculation  while  the  latter  was 
resolved  by  determining  permeability  by  the  technique  proposed 
by  Dranchuk  and  Sadiq^11^.  Following  is  a  summary  of  the  perme¬ 
ability  data. 

TABLE  2 


Core 

Liquid  Permeability 
(MDS) 

Berea  1 

4.74 

Berea  1A 

10.  33 

Berea  2 

550,4 

Limestone 

6,20 

Series 


606.0 


' 


. 


’ 
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The  scattering  of  the  steady  state  data  which  was  observed 

(23) 

at  low  flow  rates  was  also  observed  by  Mackett  in  an  earlier 
study  and  attributed  to  inaccuracies  of  the  Vol-O-Flo  meter. 
UNSTEADY  STATE  MEASUREMENTS 

A  number  of  transient  flow  tests  were  performed  on  each 
core.  The  visicorder  readings,  which  represent  pressure  dif¬ 
ferentials  across  individual  sections  of  a  core,  were  combined 
with  the  upstream  and  downstream  core-face  pressure  to  yield 
pressure  profiles  for  varying  cumulative  times. 

During  the  early  flow  period,  prior  to  the  transient 
reaching  the  external  boundary,  the  differential  pressure  readings 
were  subtracted  from  the  known  upstream  pressure  to  yield  the 
downstream  pressure  and  provide  a  check  of  that  pressure.  How¬ 
ever  on  occasion,  a  discrepancy  was  encountered  and  in  those 
incidents  the  deviation  was  distributed  proportionally  over  all 
pressure  differentials. 

At  later  times  when  the  transient  had  attained  the  upstream 
end  of  the  core  and  the  downstream  core -face  had  stabilized,  profiles 
were  obtained  by  adding  differentials  to  this  stabilized  pressure. 

Any  discrepancies  resulting  from  the  check  of  the  upstream  pressure 
were  treated  identically  to  those  encountered  at  early  flow  times. 


. 
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Negative  pressure  differentials  were  observed  during  the 
processing  of  data  for  the  Berea  2  and  Limestone  samples.  Such 
differentials  were  treated  as  zero  differentials  and  any  discrepancy 
arising  from  the  pressure  check  were  proportionally  distributed 
across  the  pressure  profile.  These  occurrences  are  explained  in 
the  following  section. 

Unsteady  state  flow  tests  on  each  sample  involving  the 
boundary  conditions  of  constant  terminal  pressure  with  constant 
pressure  at  the  external  boundary  were  continued  until  steady  state 
conditions  had  been  attained.  The  resulting  steady  state  pressure 
profile  was  utilized  to  calculate  the  permeability  of  each  section  of 
the  sample.  Table  3  lists  the  permeability  distribution  of  each 
sample.  The  permeabilities  are  in  millidarcies. 


TABLE  3 

Permeability  of  Sections 


1 

2 

Berea  1 

1.46 

42.1 

Berea  1A 

2.61 

8.05 

Berea  2 

653. 

627. 

Limestone 

10.8 

14.8 

Series 

473. 

535. 

3 

4 

5 

6 

41.2 

22.5 

14.4 

1.30 

234. 

234. 

15.0 

9.80 

674. 

655. 

476. 

378. 

283. 

282. 

2,39 

2.37 

533. 

662. 

652. 

767. 

’  . 


• 
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DISCUSSION  OF  RESULTS 

EXPERIMENTAL  RESULTS 

Experimental  unsteady  state  flow  tests  wetie  conducted  on 
the  five  core  samples  for  the  following  two  different  sets  of 
boundary  conditions. 

1.  Constant  terminal  pressure  case  with  constant 
pressure  at  the  external  boundary, 

2.  Constant  terminal  pressure  case  with  a  sealed 
external  boundary. 

A  typical  example  of  results  obtained  from  such  flow  tests 
is  presented  in  Figure  4  for  the  limestone  sample.  The  ordinate 
represents  the  pressure  term  p(x,  t),  which  is  a  function  of  dis¬ 
tance  x  and  time  t  while  the  abscissa  is  the  dimensionless  dis¬ 
tance  x/L  where  L  is  the  overall  length  of  the  system. 

The  pressure  distribution  in  Figure  4  corresponding  to 
t  -  0,1  minutes  exhibited  the  pressure  drawdown  occurring  at 
the  downstream  core-face  upon  opening  the  toggle  valve  to 
achieve  constant  terminal  pressure  at  this  end  of  the  system. 
During  this  period  the  production  originated  from  the  region  in 
the  immediate  vicinity  of  the  producing  face.  Furthermore,  as 
a  result  of  this  fluid  loss  through  production  additional  pressure 


*3 

i 
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P(x,  t),  (PSIA) 
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FIGURE  4 

PRESSURE  DISTRIBUTION 

LIMESTONE:  CONSTANT  TERMINAL  PRESSURE  WITH 
CONSTANT  PRESSURE  AT  EXTERNAL  BOUNDARY 


x/L,  (DIMENSIONLESS) 


-  29  - 


declines  occurred  more  remote  from  the  producing  end  thus  per¬ 
mitting  the  fluid  there  to  expand  and  flow  towards  this  end  of  the 
system.  In  this  way  the  pressure  disturbance  originating  at  the 
downstream  core -face  propagated  through  the  system  which  had 
behaved  as  if  it  were  infinite  as  long  as  the  disturbance  had  not 
reached  a  finite  boundary. 

Finite  behaviour  commenced  when  the  transient  reached 
the  external  boundary.  The  distribution  at  t=0,5  minutes  illustrates 
the  first  trace  of  finite  behaviour  with  constant  pressure  being 
maintained  at  the  external  boundary.  Furthermore  with  the 
commencement  of  finite  behaviour,  the  heterogeneous  nature  of 
the  core  sample  became  evident. 

The  Limestone  sample  had  a  low  permeability  section  on 
the  upstream  end  resulting  in  the  comparatively  large  pressure 
drop  observed  from  x/L  =  0.  664  to  x/L  =  1.  0.  This  tight  section 
substantially  reduced  the  flow  of  gas  to  the  more  permeable 
sections,  consequently  these  sections  depleted  almost  as  if  a 
sealed  boundary  exists  at  x/L  =  0.  664.  Upon  reaching  steady 
state  the  depletion  was  almost  complete  in  the  higher  capacity 
sections  as  exhibited  by  the  pressure  distribution. 
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The  extent  of  the  depletion  resulting  from  the  heterogeneous 
nature  of  the  sample  was  illustrated  by  comparing  the  mean  pressures 
of  the  steady  state  profiles  of  the  experimental  data  and  that  of 
a  homogeneous  system.  Had  the  core  been  homogeneous,  the 
steady  state  pressure  distribution  would  have  been  identical  to 
the  theoretical  one  depicted  by  the  dashed  line  on  Figure  4  which 
resulted  in  a  mean  pressure  of  30. 3  psia.  The  mean  pressure 
associated  with  the  experimental  steady  state  profile  was  23.6  psia. 

Application  of  the  second  set  of  boundary  conditions  to 
the  limestone  core  yielded  results  as  presented  in  Figure  5. 

Initially  the  pressure  distributions  were  identical  to  those  of  the 
first  set  of  boundary  conditions;  however,  upon  reaching  the 
external  boundary  they  were  observed  to  decline  more  rapidly. 

This  behaviour  can  be  attributed  to  the  lack  of  flow  across  the 
external  boundary.  Again,  the  tight  section  on  the  upstream 
end  of  the  sample  was  apparent. 

A  peculiar  problem  arising  while  conducting  unsteady 
state  flow  tests  on  the  Indiana  limestone  sample  was  the  occurrence 
of  negative  pressure  differentials  for  both  sets  of  boundary  conditions. 


One  possible  explanation  woul  d  be  a  localized  pressure 
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FIGURE  5 

PRESSURE  DISTRIBUTION 
LIMESTONE: CONSTANT  TERMINAL  PRESSURE  WITH 
SEALED  EXTERNAL  BOUNDARY 


x/L,  (DIMENSIONLESS) 
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expansion  in  the  immediate  vicinity  of  the  affected  pressure  taps 
resulting  in  an  decreased  fluid  velocity  and  the  resultant  negative 
differential  „ 

A  more  reasonable  explanation  might  be  a  severe 
reduction  of  the  transmissibility  of  the  pore  channels  in  the  im¬ 
mediate  vicinity  of  the  effected  pressure  taps  because  of  the 
mechanical  drilling  of  the  taps .  This  damaged  surface  would 
delay  the  response  of  the  corresponding  transducer  and  result 
in  negative  pressure  differentials , 

Assuming  the  latter  to  be  the  case,  the  affected  taps 
were  scraped  with  a  sharp  instrument  and  cleaned  out,  but  no 
improvement  was  observed.  Perhaps  if  the  core  had  been  reversed 
or  equipment  modifications  made  this  problem  might  have  been 
resolved.  Modifications  could  include  changing  to  static  pressure 
recording  transducers  or  drilling  more  than  one  pressure  tap  at 
each  length  location. 

Results  of  unsteady  flow  tests  conducted  on  the  other  core 
samples  are  graphically  illustrated  in  Appendix  D, 

Figure  6  illustrates  the  heterogeneous  behaviour  of  the 
Berea  1  core  sample.  Caehero  in  an  earlier  investigation  believed 
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FIGURE  6 

PRESSURE  DISTRIBUTION 
BEREA  1:  CONSTANT  TERMINAL  PRESSURE  WITH 
CONSTANT  PRESSURE  AT  EXTERNAL  BOUNDARY 


x/L,  (DIMENSIONLESS) 
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low  permeability  sections  on  each  end  of  this  core  were  due  to 
machining  of  the  end  faces.  In  this  study  three-eights  of  one 
inch  were  fractured  off  each  end  before  testing  and  obtaining 
results  displayed  in  Figure  6. 

A  sample  of  further  tests  performed  on  the  Berea  1 
core  but  with  additional  rock  volume  removed  from  each  end  is 
illustrated  in  Figure  7.  This  modification  more  than  doubled 
ihe  permeability  (see  Table  2)  but  did  not  entirely  remove  the 
heterogeneous  nature  of  the  core.  At  steady  state  a  large  pressure 
drop  associated  with  the  tight  section  remained  on  the  downstream  end 
while  the  upstream  end  was  considerably  improved. 

Figure  8  also  depicts  the  heterogeneous  behaviour  of 
the  Berea  1  and  Berea  1A  samples.  Earlier  assump¬ 
tions  that  this  behaviour  was  the  result  of  damage  incurred  while 
machining  the  end  face  are  disproved  by  this  investigation.  Either 
the  sample -wasnaturally  very  heterogeneous  or  it  was  damaged 
during  earlier  studies  by  a  previous  investigator. 

One  of  the  more  homogeneous  cores  tested  was  the  Berea 
2  sample.  A  typical  result  is  presented  in  Figure  9.  Although 


this  core  was  homogeneous  it  had  such  a  high  flow 
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FIGURE  7 

PRESSURE  DISTRIBUTION 

BEREA  1A:  CONSTANT  TERMINAL  PRESSURE  WITH 
CONST  ANT  PRESSURE  AT  EXTERNAL  BOUNDARY 


x/L,  (DIMENSIONLESS) 
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FIGURE  8 

PRESSURE  DISTRIBUTION 
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FIGURE  9 

PRE SSURE  DISTRIBUTION 
BEREA  2:  CONSTANT  TERMINAL  PRESSURE  WITH 
CONST  \NT  PRESSURE  AT  EXTERNAL  BOUNDARY 
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capacity,  that  constant  terminal  pressure  could  not  be  established 
prior  to  the  transients  reaching  the  upstream  end  of  the  core.  Had 
this  boundary  condition  been  achieved  instantaneously  after  opening 
the  toggle  valve,  drawdown  would  have  been  more  severe  and  a 
steady  state  condition  would  have  been  approached  more  rapidly. 
Furthermore,  upon  commencement  of  flow  at  the  external  boundary 
the  high  permeability  contributed  to  high  flow  rates  which  could 
not  be  supplied  by  the  pressure  regulator  without  a  slight  pressure 
decline.  From  an  experimental  viewpoint,  constant  terminal 
rate  would  be  a  more  appropriate  boundary  condition  to  study 
for  this  high  capacity  sample.  Such  a  boundary  condition  might 
be  achieved  by  a  Moore  constant  rate  regulator. 

Another  problem  associated  with  the  Berea  2  sample  was 
the  recording  of  negative  pressure  differentials  by  the  upstream 
transducer  while  studying  the  second  set  of  boundary  conditions. 
This  occurred  at  late  time  periods  just  prior  to  the  system  attaining 
a  steady  state  condition.  Such  a  behaviour  might  be  attributed  to 
a  system  leak  at  the  upstream  end  of  the  core  or  transducer, 
however,  neither  was  detected. 

Typical  results  obtained  while  testing  the  Series  core 
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which  comprised  of  two  different  rock  types,  an  Alundum  section  and 
a  Berea  sandstone  section  arranged  in  series,  are  presented  in  Figures 
10  and  11.  Figure  10  represents  the  unsteady  state  flow  results  when 
the  Alundun  Section  was  at  the  downstream  end  of  the  system  while 
Figure  11  are  results  for  the  same  sample  but  reversed^with  the 
Berea  section  at  the  downstream  end.  These  tests  were  intended 
to  reveal  the  effect  of  a  fracture  in  the  Berea  section  and  the  discon¬ 
tinuity  created  by  two  different  rock  types.  Unfortunately  the  taps 
were  not  spaced  closely  enough  to  detect  the  presence  of  the  fracture 
nor  were  the  permeabilities  of  the  rock  types  sufficiently  different  to 
show  the  discontinuity.  Furthermore,  Figures  10  and  11  illustrate 
that  the  high  capacity  of  the  sample  permitted  the  transients  to  reach 
the  external  boundary  before  constant  terminal  pressure  could  be 
achieved.  However,  the  tests  did  illustrate  that  the  core  depleted 
initially  much  more  rapidly  when  the  tight  section  was  at  the  down¬ 
stream  end  as  in  Figure  10  than  when  reversed  as  in  Figure  11.  At 
t=0.5  minutes  steady  state  conditions  had  not  been  attained  in  Figure  11 
whereas  they  have  in  Figure  10. 

As  a  check  of  the  accuracy  of  the  experimental  results,  a 
material  balance  based  on  the  instantaneous  pressure  profiles  observed 
was  utilized  to  yield  the  dimensionless  production  function.  The  pro¬ 
cedure  involved  numerical  integration  of  the  profiles  to  give  average 
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FIGURE  10 
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FIGURE  II 

PRESSURE  DISTRIBUTION 
SERIES  CORE:  BEREA  ALUNDUM 
CONSTANT  TERMINAL  PRESSURE  WITH  CONSTANT 
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pressures  which  when  utilized  in  a  material  balance  yielded  cumulative 

production  for  times  corresponding  to  the  profiles.  Once  obtained  these 

were  converted  to  dimensionless  production  as  defined  by  Van  Everdingen 
(31) 

and  Hurst  v  .  All  calculations  were  performed  on  the  digital  computer 
and  the  results  are  presented  in  Figure  12. 

As  the  theory  predicts,  dimensionless  production  from  all 
cores  regardless  of  its  physical  properties  are  coincidental  at  dimen¬ 
sionless  times  corresponding  to  an  infinite  behaviour  of  the  system. 
Subsequent  to  the  transients  reaching  the  finite  boundaries,  departure 
from  the  infinite  curve  seemed  to  occur  commencing  with  the  shorter 
specimens.  The  minor  discrepancies  noted  are  attributed  to  non¬ 
homogeneities  of  the  systems  studied. 

The  heterogeneous  behaviour  of  the  core  samples  is  exempli¬ 
fied  in  Figure  13.  This  type  of  dimensionless  plot  should  bring  all 
the  pressure  profiles  on  to  a  single  curve.  However,  because  there 
was  no  theoretically  valid  means  of  expressing  dimensionless  time  to 
account  for  variation  of  permeability  with  position,  deviations  from 
the  common  curve  were  observed.  The  permeability  used  in  evalu¬ 
ating  dimensionless  time  as  plotted  along  the  abscissa  in  Figure  12  and 
Figure  13  was  the  average  liquid  or  absolute  permeability  of  the  entire 
core  as  determined  from  steady  state  flow  data. 
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An  additional  check  of  the  experimental  data  was  provided 
by  a  comparison  of  the  stabilization  time  and  the  physical  properties 
of  each  system  as  presented  in  Figure  14.  The  stabilization  time 
which  was  obtained  from  visieorder  readings  represented  the  time 
that  each  and  every  point  in  the  system  had  approached  a  steady  state 
condition.  That  is,  the  time  when  the  pressure  drop  across  each 
section  had  stabilized  to  some  constant  value.  The  grouping  of 
physical  properties  as  plotted  on  the  ordinate  was  based  on  correl¬ 
ations  summarized  by  Van  Poolen  These  correlations  implied 

a  straight  line  relationship  which  seems  to  be  confirmed  by  the 
experimental  data  of  this  study.  Each  data  point  pertains  to  one 
of  the  five  flow  systems  studied. 

For  linear  systems  the  stabilization  time  can  be  predicted 
by  the  following  correlation  which  was  developed  from  this  study: 


ts 


1.31 


.344 

ft  L2 

K 


Permeability  has  the  units  of  millidarcies  and  length  has  the  units 
of  inches  while  time  is  in  minutes. 

Such  a  correlation  will  facilitate  the  planning  of  unsteady 
state  flow  tests  provided  the  physical  properties  of  the  flow  systems 


are  known. 
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FIGURE  14 

STABILIZATION  TIME 
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NUMERICAL  SOLUTIONS 

A  duplication  of  the  experimental  results  was  achieved  by 
solving  the  partial  differential  equation  describing  gas  flow  by  the 
numerical  technique  presented  in  an  earlier  section  and  developed 
in  the  Appendices, 

Accuracy  of  the  numerical  solution  for  a  homogeneous 
system  is  illustrated  in  Figure  15  through  a  comparison  with  an 
analytical  solution  presented  in  Katz*18\  Discrepancies  between 
the  numerical  and  analytical  solutions  show  a  trend  of  going  from 
positive  deviations  to  increasingly  more  negative  deviations  as 
time  increases,  A  further  comparison  could  not  be  made  because 
the  pressure  profiles  have  reached  the  external  boundary  and  the 
analytical  solution  applies  only  to  an  infinite  system.  A  comparison 
performed  on  data  for  the  other  samples  revealed  similar  discrep¬ 
ancies;  however,  there  was  no  consistency  as  to  whether  the 
numerical  solutions  were  greater  than  or  less  than  the  analytical 
solution  for  a  corresponding  time.  If  the  discrepancies  were 
consistently  in  one  direction  they  might  suggest  some  inherent 
error  in  the  mathematical  model. 

Utilizing  the  average  Klinkenberg  permeability  of  each  core 
and  the  mathematical  model,  pressure  profiles  for  homogeneous 
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FIGURE  15 
BEREA  1 
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systems  were  calculated.  One  such  profile  is  illustrated  in 
Figure  16.  The  ordinate  is  a  dimensionless  pressure  term 
P(xjf)/Pj?  which  is  the  ratio  of  pressure  change  at  any  point 
within  the  porous  medium  to  the  constant  pressure  maintained  at 
the  external  boundary.  The  abscissa  is  the  dimensionless  distance, 
x/L,  Unlike  the  experimental  data  presented  in  Figure  13,  a 
single  curve  was  obtained  with  the  exception  of  that  for  the  Berea 
2  core.  Figure  16  supports  the  theoretical  aspects  of  such  a  plot 
and  illustrates  the  effects  of  heterogeneous  systems  when  compared 
with  Figure  13. 

The  Berea  2  data  deviates  from  the  common  curve  because 
the  constant  terminal  pressure  boundary  condition,  common  to 
the  other  tests  was  not  attained  immediately,  consequently  at  a 
corresponding  dimensionless  time,  depletion  was  less  and  the 
pressure  profile  was  resultantly  higher. 

Extending  the  numerical  technique  to  the  non-homogeneous  systems 
by  using  the  actual  permeability  of  each  section  studied  experimentally, 
results  in  the  comparison  presented  in  Figure  17.  The  numerical 
approximation  coincided  with  the  experimental  results  for  early 
times  but  deviates  positively  in  the  interior  sections  as  steady 
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FIGURE  16 
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FIGURE  17 

PRESSURE  DISTRIBUTION 

BEREA  1A:  CONSTANT  TERMINAL  PRESSURE  WITH 
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state  conditions  were  approached.  At  flow  periods  as  short  as 

t=0.2  minutes,  it  was  apparent  that  the  depletion  of  the  system  as 

predicted  by  the  numerical  solution  was  less  than  that  experimentally 

(10) 

observed.  Chwyl  who  by  a  numerical  technique  was  able  to 
examine  the  contributions  of  gas  slippage  and  inertia  individually 
also  observed  this  phenomena.  Molecular  streaming  or  gas 
slippage,  which  was  not  accounted  for  in  the  mathematical  model, 
was  the  most  apparent  cause  of  such  behaviour.  However  utilizing 
Chwyl  fs  data,  slippage  accounted  for  only  about  0.1  psi  of  the 
discrepancy.  Using  the  same  source,  inertial  effects  were  calcu¬ 
lated  to  be  negligible. 

Similar  deviations  are  observed  when  comparing  results 
for  the  second  set  of  boundary  conditions  as  applied  to  the  Berea 
1A  sample  and  presented  in  Figure  18. 

Figure  19  compares  results  for  the  Berea  2  sample  which 
experienced  the  highest  flow  rates  encountered.  For  corresponding 
time,  the  experimental  data  were  greater  than  the  numerical  solution, 
just  opposite  to  the  previous  comparison.  This  behaviour  can  be 
explained  by  the  fact  that  higher  internal  pressures  were  necessary 
to  provide  sufficient  pressure  drop  to  overcome  the  inertial  effects. 
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FIGURE  18 

PRESSURE  DISTRIBUTION 

BEREA  i  A:  CONSTANT  TERMINAL  PRESSURE  WITH 
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FIG  UP  E  19 
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not  considered  in  the  mathematical  model.  Evaluation  of  the 
dimensionless  inertial  resistance  coefficient  and  the  dimensionless 
slip  coefficient  as  defined  by  Chwyl,  revealed  these  inertial  effects 
might  be  significant.  However,  their  actual  quantitative  contribu¬ 
tion  towards  the  discrepancies  observed  in  Figure  19  could  not  be 
evaluated  because  of  the  exponential  behaviour  displayed  by  Chwyl’s 
results  and  the  need  to  interpolate.  Nevertheless,  from  qualitative 
aspects  it  appeared  that  ChwyPs  results  would  not  account  for  the 
entire  disrepancy  observed. 

Additional  factors  which  would  contribute  to  the  discrepancies 
observed  in  Figures  17,  18  and  19  are  erroneous  permeabilities 
and/or  an  incorrect  mathematical  model.  Possible  errors  in  the 
permeability  data  arise  because  of  the  limited  amount  of  steady 
state  data  available  to  evaluate  the  variation  of  permeability  with 
position.  In  the  case  of  the  mathematical  model,  the  method  of 
discretizing  the  derivatives  resulted  in  an  equation  which  produced 
a  non-symmetrical  matrix  of  non-linear  terms.  Such  a  matrix  could 
yield  cyclic  results  which  may  be  the  situation  in  Figure  18. 

Experimental  and  numerical  solutions  for  the  other  core 
samples  are  presented  in  Appendix  D,  Figures  D-8  through  to 
Figure  D-10. 
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CONCLUSIONS 

From  observations  made  and  results  obtained  in  the  course 

of  this  investigation,  the  following  conclusions  were  formulated, 

1.  The  mathematical  model  developed  for  linear  systems  can 
be  solved  numerically  to  give  results  approximating  those 
experimentally  observed.  Deviations  between  results  can 
be  attributed  to  molecular  streaming,  inertial  effects, 
erroneous  permeabilities  and/or  inherent  errors  in  the 
model. 

2.  No  valid  means  of  expressing  dimensionless  time  to 
account  for  variation  of  permeability  with  position  in 
a  linear  system  exists. 

3.  Pressure  distributions  obtained  experimentally  for  the 
two  sets  of  boundary  conditions  are  identical  for  each 
system  until  the  pressure  disturbance  created  at  the  down¬ 
stream  end  reaches  the  upstream  end. 

4.  The  Series  sample  designed  to  simulate  non -homogeneous 
behaviour,  is  not  heterogeneous  enough  and  possesses  too 
high  a  flow  capacity  for  optimum  transient  study, 

5.  The  excessive  pressure  differentials  observed  across 
each  end  of  the  Berea  1  sample  are  not  attributable  to 
mechanical  end  effects  but  arise  as  a  result  of  damage 
incurred  in  a  previous  investigation. 
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Long  core  samples  are  desirable.  Lengths  in  excess  of 
five  feet  facilitate  the  measurement  of  transients . 

Core  samples  with  an  absolute  permeability  in  excess  of 
100  millidarcies  are  undesirable  for  transient  flow  tests. 
This  type  of  sample  makes  it  difficult  to  simulate  the 
boundary  conditions  examined  in  this  study. 


RECOMMENDATIONS 


More  representative  transient  pressure  measurement  might 
be  obtainable  by  drilling  more  than  one  pressure  tap  at  each 
location  along  the  length  of  the  flow  system. 

Numerous  steady  state  flow  tests  should  be  performed  on 
each  core  sample  to  facilitate  the  evaluation  of  the  per¬ 
meability  of  each  section  of  the  sample. 

When  estimating  time  required  to  conduct  transient  flow 
tests,  the  correlation  developed  in  this  study  relating 
stabilization  time  to  physical  properties  should  be  employed. 
The  boundary  condition  of  constant  terminal  rate  should  be 
investigated  because  it  has  more  practical  applications  and 
it  might  facilitate  the  transient  study  of  higher  capacity 
samples . 

The  mathematical  model  should  be  modified  to  incorporate 
molecular  streaming  and  inertial  effects,  and  perhaps  obtain 
a  better  approximation  to  the  experimental  results. 
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NOMENCLATURE 

[a]  matrix  of  non-linear  terms 

b  vector  of  known  terms 

K  absolute  (liquid)  permeability  (L2) 

2 

k  permeability  (L  ) 

kg  gas  permeability  (L  ) 

L  overall  length  of  linear  system  (L) 

M  molecular  weight  (m) 

p  pressure  (m/Lt2) 

p  vector  of  unknown  pressures  (m/lt2) 

Pj  initial  pressure  (m/Lt  ) 

Pw  pressure  at  the  producing  end  (m/Lt2) 

R  universal  gas  constant  (mL2/t2T) 

t  time  (t) 

ts  stabilization  time  (t) 

T  temperature  (T) 

v  volumetric  flow  rate  per  unit  cross-sectional  area  (L/T) 

W  inverse  of  product  of  viscosity  and  compressibility  (Lt ) 

m 

x  distance  along  linear  axis  (L) 

z  gas  compressibility  factor 

n  viscosity  (m/Lt) 

jLig  gas  viscosity  (m/Lt) 
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density  (m/L^) 

Q 

gas  density  (m/L  ) 
porosity 
SUBSCRIPTS 
g  gas 

i  initial  value  or  conditions 

j  grid  point  index  along  length 

k  grid  point  index  along  time 

s  stabilization 

w  producing  end  of  the  system 

SUPERSCRIPTS 
-*■  denotes  vector  quantity 

r  index  of  iteration  number 

denotes  a  mean  value 
MATHEMATICAL  NOTATION 

7  del  operator  (gradient) 

A  delta  difference 

V*  divergence 

y  partial  derivative  with  respect  to  y 
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DERIVATION  OF  GAS  FLOW  MATHEMATIC AL  MODEL 

The  continuity  equation  for  a  compressible  fluid  may  be 
written  as: 


Amount  of 

Amount  of 

Increase  of  Mass 

Mass  Input 

Mass  Output 

Within  the  Region 

or  in  mathematical  terms,  this  material  balance  based  on  conser¬ 
vation  of  mass  for  isothermal  fluid  flow  through  a  porous  medium 


is: 


v.[>>  - 

3 1 


(A-l) 


Assuming  flow  to  be  laminar  and  horizontal,  the  velocity 
vector  in  Equation  (A-l)  -may  be  expressed  by  Darcy’s  law: 


v=  ~k(P)  \7P 
*(P) 


Substituting  this  into  Equation  (A-l)  for  the  velocity  vector 

yields: 


(A-2) 


Since  the fiLuid  being  considered  is  a  real  gas,  the  density  terms 

in  Equation  (A-2)  may  be  represented  respectively  by  the  following 

\ 

forms  of  the  equation  of  state: 
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and 


dPg_  =  m_ 
Bt  RT 


3  P 
2)  t  a(p) 


which  when  substituted  into  Equation  (A-2)  gives: 


V- 

-M 

(kS(p)p  Inp] 

ii 

i 

£ 

cv 

p 

RT 

(ng(p)  s(p)  J  _ 

ht  at 

_*■  (p)  _ 

Since  M/RT  is  constant  for  isothermal  variation  of  pressure  it  can¬ 
cels  as  do  the  negative  signs  yielding: 


V* 

kg(P)  -i 

- 5 - pVp 

=  jzf  d 

p 

big(p)  g(p) 

3t 

_z(p)_ 

(A-4) 


In  this  study  involving  very  heterogeneous  cores,  the  perme¬ 
ability  was  considered  a  function  of  position  only,  independent  of  any 
pressure  variation.  Furthermore,  the  systems  were  linear  therefore 
the  resulting  equation  is: 


K(x)  p 

~d  P 

=  $  3 

p 

_^g(p)  B(p) 

"d  x 

dt 

_*(P)  _ 

If  a  pressure  function  is  introduced  such  that: 

W(p)  =  1 

JUg(p)  2(p) 

Equation  (A-5)  reduces  to  a  non-linear  partial  differential  equation 


of  the  form: 


■ 
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(A-f) 


Further  modificatjcns  involve  the  introduction  of  a  mean  pressure 
term  on  the  right-hand  side  of  Equation  (A-6)  and  utilizing  the 
equivalents 


and  2p  9  p  -  d 

' .  :  dt  at 


3  x  9k 


to  yield: 


3_  W(p)  K(x)3p2  =  J)T p2 


(A-7) 


^ x  —  9  x_  p  9t  _g(p) 


An  appropriate  set  of  boundary  and  initial  conditions  are: 
1.  Constant  terminal  fires  sure  case  with  external  boundary 

i* 

pressure  constant 


INITIAL  CONDITION^ 
p(x,  o)  =  pi 

BOUNDARY  CONDITIONS 

P  (o,  t)  =  Pw; 

* 

P  CM)  =  Pi 


2 .  Constant  terminal  pressure  case  with  external  boundary 


sealed 
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INITIAL  CONDITION 

p(x,o)  =  p. 

BOUNDARY  CONDITIONS 
P(°,t)  =  Pw 


d  p  (L,t)  =  o 
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DISCRETIZATION  OF  PARTIAL 
DIFFERENTIAL  EQUATION  AND  METHOD  OF  SOLUTION 

CRANK-NICHOLSON  IMP  LICIT -ONE  DIMENSION 


From  Appendix  A,  the  diffusivity  equation  for  flow  of  a  real 


gas  is: 


3 

W(p)  K(x)  <?p2 

II 

CD 

2~| 

P 

dx 

P  ^t 

_z(p]_ 

Differentiating  Equation  (B-l)  gives: 

K(x)  <?p2  <?W{p)  +  W(P)  K(x)  a2p2  +  d p2 W(p)  i?K(x) 
dx  dx  5x2  dx  dx 


(B-l) 


d 

p  77 


p 


_z(P)_. 


(B-2) 


To  formulate  a  difference  equation  for  Equation  (B-2)  a  network 
having  mesh  width  of  Ax  and  At  is  established  as  portrayed  in  Figure 


(B-l). 


TaylorTs  series  expansion  provides  the  approximation  to 


the  spatial  derivatives.  Consider  any  point  j  on  the  grid  spacing 


P1+i  y  SP,  L  +  dv*  Ax  +  dZp2  Ax^  +  .  .  ..  (B-3) 
3  1,k  Lk  ±  —^3—  ~2l 


Pj-l,k  =  Pj,k--^E-Z  X+%J? 

ox  o  x*  21 


j(B-4) 


Subtraction  of  Equation  (B-3)  and  Equation  (B-4)  neglecting  terms 

3  3 

of  the  order  of  d  p  Ax  and  higher  gives: 

7?  31 


•  ;  i.  h 
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2  2 
P  j+l,k  -  P 


k 


2  Ax 


(B-5) 


The  addition  of  Equation  (B~3)  and  Equation  (B-4)  with  terms  of 


4  2 

order  a  p 


A  x  and  higher  neglected  yields : 


x 


4! 


q2  2 

d  x 


P  j+l,k  2pj,k+pj-l,k 
~Z - - 


Ax 


(B-6) 


Similarly  expanding  p.  , 

j,K+I 


about  p,  , 
j.k 


and  neglecting  the  second  and 


higher  order  terms  gives: 


r  2n 
p 

=  i 

pj,k+l  _  Pi.k  “1 

b(p)J 

At 

LZ(Pj>k+1)  z<pj>k>J 

(B-7) 


This  derivative  will  be  centered  at  k-hl/2. 

Finally,  the  derivative  dK(x}/  d  x  differs  from  the  derivative 
of  Equation  (B-5)  only  in  the  respect  that  values  of  the  property 
K(x)  are  located  at  Therefore  the  derivative  approximation 

y 

is 

=  A^L  =  K<xt-rt/2)  -  (B-8) 

Equations  (B-5)  and  (B-6)  represent  the  finite  difference 


approximations  to  the  spatial,  derivatives. 


•  .  •  ’  . 
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However  since  it  is  desirable  to  employ  Crank-Nicholson  approxi¬ 
mations,  the  appropriate  equations  are: 


3  2 

O  P  = 


4  Ax 


P 


2  2  2  — I 

j^,K"  p  j_i,k  +P  j+l»k+l  ~p  j-l,k+l| 


(B-9) 


«2  2 

1_P2  =  __1 

2  Ax2 


P2j+l.k  -  2P2J,k  +p2]-V  +PVl»k+l  '  2p2j,k+l  +p2j- 

(B-10) 


The  substitution  of  Equations  (B-7),  (B-8),  B-9)  and  (B-10)  into 
Equation  (B-2),  and  the  collection  of  terms  according  to  their  time 
and  space  subscripts  with  respect  to  pressure  yields: 

|k(Xj)  AWtej)  +  2W(pj)  AK(xj)  +4W  fty  K(Xj)l  p2j+i 


+  8^  Ax2 


z(pj)  “  11.  •  L 


2 

p  j+i,k 


-  SW^)  K(Xj|  P  j(k  +  |  K(Xj)  AW(Pj)  -  2W(pj)  AK(Xj) 

2 


■KCxp  AW(p  )  -  2W(Pj)  AK(Xj)  -  4W(pj)  K(Xj) 


j+l*k-d 


+ 

2 

p 


8$  Ax2 


P  At  x(p 


jJ 


+  8W(p  )  K(x  ) 


!j,k+i  +  p(xJ>  4W(pj)  +2W(pj)  4K<x 


j-l,k+l 


(B-ll) 


,  i 


l,k+l 


+  4W(pj)K(x„) 


.)  -  4W(p3)K(x^) 
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wherein 

AK(Xj)  =  K(xj+1/2)  -  K(^_1/2) 

AW(pj)=  W(pj4)  _  W(Pj_j) 

SOLUTION  OF  THE  NON-LINEAR  C-N  IMPLICIT  EQUATION 

The  C-N  implicit  procedure  in  one-dimension  as  illustrated 
in  Figure  (B-l)  was  applied  utilizing  Equation  B-ll  with  appropriate 
initial  and  boundary  conditions  to  obtain  the  pressure  distribution 
for  each  time  step. 

For  convenience  Equation  (B-ll)  is  written  as  : 

A^j-l,k^)  +B<P2j)k+1>  +  C<P2j+l,k*>  -  Dj(k  (B-12) 

wherein: 

A=  K(xj)  AW(p.)  +  2W(p .)  AK(x.)  -  4W(p .)  K(Xj) 

B=  8^ Ax2  +  8W(p  )  K(x.) 

pAt  z(pj)  3 

C=  -K(x.)AW(p.)  -  2W(Pj)AK(x.)  -  4W(p.)  K(Xj) 


and 


' 
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Time  (t)  Sweep  Direction  - ► 


j  1  2  3  4  5  n-  2  n-1  n 

x=0  x=L 

ONE  DIMENSIONAL  LINEAR  CRANK  NICHOLSON  IMPLICIT 
PROCEDURE  GRID  POINTS  AND  TIME  LEVELS 


FIGURE  B-l 
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D= 


K(xj)AW(pj)  +2W(Pj)AK(x.)  +  4W(p.)  K(Xj) 


P  H,k 


4~MAx2  -  8W(p  )K(x.) 

A  -f  /«  V  J  J 


jp  Atz(pj) 


2 

P  j,k 


-K(Xj)AW(p.)  -  2  W(p.)AK(x.)  +4W(p,)K(x.) 


P  j-l.k 


Solution  of  the  pressure  distribution  at  time  step  k+L 
commenced  by  guessing  a  pressure  profile  at  this  time  step  to 
permit  the  evaluation  of  A,  B  and  C.  The  steady  state  equation 
provided  this  initial  guessed  profile.  Subsequent  distributions 
employed  to  evaluate  A,  B  and  C  represent  solutions  obtained  in 
a  previous  iteration  step. 

The  application  of  Equation  (B-12)  to  the  first  set  of 
boundary  conditions  is  as  follows: 

GRID  POINT  (j=2)  —  FIGURE  (B-l-a) 

At  this  location  Equation  (B-12)  takes  on  the  form: 

(A)  R2  M  +  (B)  p2  k+1  +(C)  p23i  k+1  =  Dj  k  CB-18) 


however  from  boxmdary  conditions  presented  in  Appendix  A 

constant  terminal  pressure  implies 

2  2 
Pi,  1^4-1  Pw 


■ 
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therefore  Equation  (B-13)  becomes 

(B)  P2,k+1+<C>  P23>k+1-  Dl,k-(A>Pw  <B'14> 

2  2 

where  F^,k+1  an<^  ft)  k-ri  are  on^y  unknowns,  since  B  and  C 
were  evaluated  from  the  guessed  profile . 

INTERIOR  GRID  POINTS  (2<j<n-l)  —  FIGURE  (B-l-b) 


For  all  interior  grid  points  Equation  (B-12)  becomes: 

(A)  P2H,k4  +  (B)P2j,k+l  +(C)  P2j+1,k+l  -  DJ  k  (B-15) 


2  2  2 

where  p  j+l  k+1  ,  P  j  and  p  ._l  k+l  are  unknowns 

GRID  POINT  (j=n-l)  —  FIGURE  (B-l-«) 

The  external  boundary  condition  states 

2  _  „2 
p  n,k+l~p  i 

therefore  Equation  (B-12)  at  this  location  is: 

(A)  pV2, k4  +  (B)  P2n.l,  k4  =  (D)n-l,  k  -  (C)p2i  (B-16) 


The  one  relationship  of  Equation  (B-14),  the  n-2  relationships 
of  Equation  (B-15)  and  the  one  relationship  of  Equation  (B-16)  provide 
for  n  conditions  for  evaluating  improved  values  of  p.  k+p  2^n-l, 
over  t^bfee  guessed.  For  example,  if  n=6,  these  conditions 


...  ■ 
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are  as  represented  in  Table  B-l 


TABLE  B-l 

C-N  IMPLICIT  PROCEDURE 


j  PX  P2  P3  P4 


p5 


p6 


1 

2 

3 

4 

5 

6 


a2 


c2 

b3 

a4 


C5 

b6 


or  in  matrix  notation 

[a]  ?  -  b 


2 

D1  "  Ai  Pw 
D2 

D3 

D4 

n5  n  2 
D6  ”  C6  Pi 


(B-l  7) 


Solutions  of  the  pressure  vector  p  was  obtained  by  matrix 
inversion  utilizing  the  square -root  algorithm  modified  for  non- 
symmetrie  matrices.  Hence 


P 


=  [a"3  t 


Evaluating  the  second  set  of  boundary  conditions  involving 
the  sealed  boundary  required  altering  only  two  of  the  n  equations. 
The  boundary  condition  implies : 


^  P  n9  k+1  -  o 

9  x 

which  can  be  incorporated  into  Equation  (B-13)  by  setting 

2  2 
P  n+lsk+l  ~  P  n-l,k+l 


(B-21) 


! 


* 
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Interior  grid  points,  to  which  Equation  (B-15)  applies  now  extend 
over  the  interval  2<j<n  and  at: 

GRID  POINT  (j  =  n)  —  FIGURE  (B-l-d) 

Introducing  the  external  boundary  condition  as  stated 
by  Equation  (B-21)  to  Equation  (B-13)  yields 

<A^>  P2  n-l,  k+l  +  «*>  P'n, k+l  =  <D>  n, k  (B"22> 

From  this  stage  on  the  solution  of  the  pressure  profiles 
is  identical  to  that  employed  in  solving  the  first  set  of  boundary 
conditions . 

Because  the  partial  differential  equation  describing  gas 
flow  is  non-linear,  it  is  necessary  to  iterate  at  each  time  step. 

If  r  indicates  the  rth  iteration,  a  solution  for  p  was  considered 
to  be  obtained  when: 


Where  £  represents  the  maximum  permissive  difference  between 

pressures  of  successive  iterations.  In  this  investigations  was 

-3 

set  equal  to  2  x  10  so  that  the  difference  between  the  accepted 
pressure  and  the  pressure  at  the  previous  iteration  was  less 
than  0.1  pounds  per  square  inch. 


. 
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COMPUTER  PROGRAM 


The  computer  program,  which  was  written  to  solve  the 
unsteady  state  gas  flow  problem  for  the  set  of  boundary  conditions 
of  constant  terminal  pressure  with  constant  pressure  at  the  external 
boundary,  may  be  divided  into  six  parts. 

(1)  MAIN  PROGRAM  -  The  main  program  computes  the 
pressure  distribution  in  a  linear  heterogeneous  system 
by  solving  the  continuity-Darcy  flow  equations  for  gas. 

A  Crank-Nicholson  implicit  procedure  is  employed  along 
with  a  Function  iterative  technique. 

(2)  PRESSURE  COEFFICIENT  SUBROUTINE  -  The  elements 
of  the  matrix  at  time  k-41  and  the  elements  of  the  vector 
corresponding  to  the  known  pressure  distribution  are 
calculated  by  this  subroutine  as  requested  by  the  main 
program. 

(3)  COMPRESSIBILITY  SUBROUTINE  -  Utilizing  Lagrangian 

/- 1 

iterpolation  and  correlations  of  Helsenrath  et  al  ,  this 
subroutine  provides  the  subroutine  discussed  in  Part  (2) 
with  compressibility  of  nitrogen  gas  as  a  function  of 
pressure  and  temperature. 

(4)  VISCOSITY  SUBROUTINE  -  This  subroutine  employs  a 
correlation  proposed  by  Kestin  et  al  to  calculate  the 


- 
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viscosity  of  nitrogen  gas  for  SUBROUTINE  PRCOEFF, 

(5)  EFFECTIVE  PERMEABILITY  SUBROUTINE  -  The  varia 
tion  of  permeability  with  distance  is  provided  by  this  sub 
routine  for  the  subroutine  discussed  in  Part  (2). 

(6)  SQUARE-ROOT  MATRIX  INVERSION  SUBROUTINE  - 
This  subroutine  inverts  the  matrix  developed  by  the  sub¬ 
routine  whose  function  was  described  in  Part  (2).  In 
matrix  notation  the  equation  describing  unsteady  state 
gas  flow  at  every  grid  point  takes  on  the  following  form, 

[a]p=“B 

Firstly  the  subroutine  inverts  the  matrix  [Aj, 

and  then  solves  for  the  vector  |T 
or  [l]^=[Ar1b- 

p=t¥1b 

which  represents  the  unknown  pressure  distribution  at 
time  k+1. 

To  solve  the  second  set  of  boundary  conditions  the  main 
program  must  be  modified  by  altering  the  equation  which 
evaluates  the  vector  b.  The  only  contribution  to  this 


Vb 
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vector  from  the  boundary  conditions  will  be  from  the 
constant  terminal  pressure  boundary.  Furthermore,  to 
SUBROUTINE  PRCOEF  must  be  added  an  equation  which . 
evaluates  the  element  at  the  (n+l)Ax  position.  These 
alterations  are  also  explained  in  Appendix  B. 
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COMPUTER  PROGRAM  NOMENCLATURE 

The  following  list  contains  those  quantities  which  are 
pertinent  to  the  understanding  of  the  program. 


A 

matrix  of  non-linear  terms  at  (i+l)^t 

A 

A 

elements  of  jAj  at  the  boundaries 

ALIP 

Lagrangian  interpolation  polynomial 

B) 

C) 

D) 

!  v  > 

coefficients  of  equation  describing  compressibility 
factor 

BA 

vector  of  known  quantities  at  time  (i)At 

BBB 

slip  coefficient 

BC 

matrix  at  time  (i)At 

BELT 

time  increment 

DELTEC 

temperature  differential  for  viscosity  correction 

DEL  WPS 

differential  of  WPS  from  one  grid  point  to  the  next 

DELXSQ 

distance  increment  squared 

DIFP 

factional  difference  between  pressures  for  two 

sequential  iterations 

DP 

DPRESS 

dime  nsi  onal  p  res  sure 

DT 

dimensionless  time 

DX 

dimensionless  distance 
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EFFK 

IA 

IG 

n 

JA 
JJ 
LA 
LK  ) 

I  ) 

LL 

P  ) 

PA) 

PERM 

PERML 

PHI 

P4  ) 
PJ|  ) 

PJPI1) 

PJPI£) 

PJPI3) 

PM 

PN  ) 
PX  ) 


effective  permeability  at  grid  point 
number  of  time  steps 
number  of  sets  of  B,  C,  D 
permeability  flag 

number  of  grid  points  at  any  time  interval 
matrix  flag 
number  of  data  sets 
time  step 

number  of  rows  and  columns  in  matrix  [Af 
pressure 

liquid  permeability  of  distance  increment 
liquid  permeability  of  core  sample 
porosity 

components  of  elements  making  up  matrix 


mean  pressure 

iterated  pressure  evaluated  for  time  (i+l)At 


• 

. 

• 
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PNSQ 

PN  squared 

PO 

known  pressure  including  initial  conditions 

T 

time 

TE 

temperature 

TP 

temperatures  associated  with  B,  C,  D 

URATIO 

ratio  of  gas  viscosities  at  P  and  one  atmosphere 

vn£ 

viscosity 

WPRS 

inverse  of  product  of  compressibility  and  viscosity 

X 

distance 

Z 

compressibilities  corresponding  to  TP 

ZN2 

interpolated  compressibility  factor 

The  preceding  nomenclature  summary  does  not  include 
that  used  in  subroutine  SQROOT .  This  subroutine  was  modified 
and  utilized  by  Allada,  and  only  the  following  must  be  known  to 


use  it: 

A 

the  matrix  of  nonlinear  terms  |AJ 

BA 

the  known  vector  b 

M=LL 

the  size  of  the  matrix  [A] 

•Y 

the  solution  vector  PNSQ 

jflibfloqas-noo  me 


'’I  3  .  i:  a  t  dr  j'f) 
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c 

c 

c 

CONSTANT  TERMINAL  PRESSURE  CASE 

WITH  CONSTANT  PRESSURE  AT  THE  EXTERNAL  BOUNDARY 

DIMENSION  A (20, 20)  *BA(  13)  ,80(20,20)  ,  DPRE  SS  (  20 , 20  )  ,  PN(  1 

IS) ,PNSQ( 13) , P0( 20, 15)  ,PX< IS ) ,DT( 15) , DPI  15, 15) »DX< 15 ) 
COMMON  HI  10)  .C(10)*D(L0}*P(15),PM(15)»T{15)  »  X  (  1 5  )  »  T  P  (  1 

c 

10) , PERM (13) ,  BBB 

COMMON  IG,  JJ » JA fLL , LK 

c 

200 

READ  IN  NUMBER  OF  DATA  SETS 

READ!  5,20C)LA 

FORMAT ( iX  ,  I  3  ) 

c 

c 

READ  IN  TEMPERATURE  VALUES 

R  E  A  D  (  5  , 106) IGf ITPI I ) , I  =  1 »  I G  ) 

c 

106 

FORMAT (lX»i3»6F6.2 ) 

c 

READ  IN  COMPRESS  I BL I T Y  CONSTANTS 

1 G  1 

READ  (5 ,101  )  (BiI),C(n,0(I),I  =  l,6) 

FORMAT ( 1X,6E10.4J 

DO  6  0  M  = 1 1  L A  . 

c 

201 

WRI TE( 6 , 201 ) M 

FORMAT!  MCCRE  N0.=',I3) 

c 

LOO 

READ  IN  BASIC  ROCK  PROPERTIES  AND  CQNDT IONS 

RE AD ( 5 , 100)  PHI,PERML,TE, IA, JA 

F  G  a  a  .A  I  (1X.F5^3,3X,F6 .2  ±  3  X  ,  F  5 . 1,2(3  Xt.I  3 1) _ _ 

115 

WRITE (6,115) PHI , PER ML , TE,  1A,  JA 

FORMAT  l  *  1  ,F5.3,3X,F6.2,3X,F5. 1,2( 3X, 13) ) 

JAM1=JA-1 

107 

RE AD (5, 107)  (PERM1J)  ,J  =  1,  JAM1) 

FORMAT! IX, 10F7.2) 

M  !ITE<6, 10  7)  (PER  M  (  J  )  ,  J=1,JAM1_L. 

108 

READ ( 5,  I  08 ) BBB 

FORMAT  (  IX,)- 7. 3) 

WRI TE ( 6 ,  108 ) BBB 

c 

c 

READ  IN  bOUNDARY  CONDITIONS 

READ ( 5  » 102 ) <  PG ( 1 >  J )  ,J=1,JA) 

102 

FORMAT ( IX* 7F10. 3) 

WRITE (6, 102)  (P0(1, J)  ,  J=l, JA) 

READ! 5,  103) ( P0(  I , 1 )  , PO ( I , JA )  , 1=2 , I  A ) 

c 

103 

FORMAT l 1X,7F10.3J 

WRITE (6, 103) ( P0( 1 , 1) ,P0( i ,JA) , 1  =  2, I  A) 

c 

104 

READ  IN  TIME  VALUES 

RE AD ( 5 , 104) (T(I)  ,  1  =  i  »  I A  ) 

FORMAT ( IX, ILF7.3) 

WRITE ( 6, 104) ( T( I) » 1=1 f IA) 

L 

c 

READ  IN  DISTANCE  VALUES  . . 

105 

READ! 5 , 105) (  X  (  J  1  *  J=  1 » JA) 

FORMAT (  LX  ,  13F6.2  ) 

WRITE(6, 105) ( X( J ) , J=1 , JA) 

c  _  1  J  V  :  t  T  IK.*.  I  fcHJJ 

I  J  i  c  I  fc  J 
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1AM 1=1  A— 1 
DO  40  I=1,1AM1 

.  LK=  I . . . . . . . . .  . 

C 

C  CALCULATE  COEFFICIENTS  OF  PRESSURES  AT  TIME=T+1 

_ DO  5  J=1,JA _ 

5  PN(  JJ^SQRTi  l  XI  JI/XI  JA)  )*{  PO  U  +  l,  J  A)*PCJ(  I  +  l,  JAJ-PCH  1  +  1, 

m^po  d  +  i,in  +pc(  i  +  i,i  )*pg(  i+i,in 

. 6  DO. . Z. «J..— JL..t  U. A. - _ - - — - - - ~ - 

P  (  J  )  =  P  N  (  J  ) 

7  PM( J)=.5*(PN( JJ+POI It J) ) 

_ JJ=1 _ : _ 

CALL  PRCO EF ( PH  I , T  E , A ) 

c 

JL _  CALCULATE  COEFFICIENTS  OF  PRESSU RES..... 41 . I IAE-5.L-. . . . . . 

DO  10  J=1,Ja 
10  P(J)=PO(I,J) 

_ AAz2 _ 

call  prcoemphi ,te,bo 
c 

JCL _ _ CA1XJJLLA.I  L.  _ _ _ _ _ 

DO  15  L=2,JAM1 
BA ( L ) =0 . 0 
DO  15  J—  1  ,  J  A 

A  1  =  0  •  0 

IF ( L  .EQ.2  .AND. J  .EU. 1 ) A1  =  A( 2 , 1 ) 

_ _ _ _ _ _ 42  =Q  ♦  0 - - - - - - ...... . . . . . . . . . . . . . . 

IF(L.EQ.JAM1.AND.J.£Q.JA)A2=A(JAM1,JA) 

15  BA(D  =BA(  L)  +  BC(L,  J)*P(  J  )*?i  J)-A1*P0(  I +1  ,  1)*P0(  1  +  1,  1 )  —  A 

_ 1 2*P0 ( 1  +  1, JA)*PO( 1  +  1, JA) _ 

J  A  M  2  =  J  A—  2 
LL= J AM2 

_ _ _ _ DO  16  L=  1 ,  JA  M2  _ 

DO  16  J= 1 »  J ALL 

16  A ( L , J) =A( L+l , J+i ) 

_ DO  17  L=1.JAM2 _ _ 

17  B A { L ) =BA ( L  +  I  ) 

C 

C  ,  CALCULATE  PRESSURES  AT  TIME=T+1 . _ . _ . . . .... 

CALL  SQRGGT ( A , B A , PNSQ ) 

C 

_C _ CHECK  CONVERGENCE  OF  PRESSURE  VALUES _ 

I  1  =  0 

DO  20  J  =  2  ,  J AMI 

_ _ PX(  J)=SQRT(PNSQiJ-m . _ . . 

0  I  FP  =  ABS ( (PX(J)-PN(J)  )  / PX ( J ) ) 

IFIDIFP.GI .0.002)11=1 

_ PN ( J )  =  PX ( J  ) _ 

20  PO ( I  +  1 , J ) =PN ( J  ) 

I  F  (  II.EQ.DGO  TO  6 

40  CONTINUE  _ ___ _ _ 

WR I T  E ( 6 , 1 1 0 ) 

110  FORMAT ( '-PRESSURE  PROFIlE  FOR  CONSTANT  PRESSURE  AT  E», 

internal  bcunlakym 
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UO  45  1=1, I A 

45  WRlTEIo, 111)T{ I ) , (PO( l , J) , J=l, JA) 

ill  F  0  R  M  ATC  « >F7.3.2X. 13 (1X.F6  *21 J _ _ _ 

DO  5  7  1=1,  I A 
DO  57  J=1,JA 

1F(  (PU(  I  *  1>**2. )~i PQ(  I  , JA)**2. ) ) 55,56,55 _ 

55  DPRESSU  ,  J  i  ==  C  PO  (  I ,  J)**2.-P0(  1  ,  JA)**2.  )/ < PO( I , 1) **2 .-PO 
1  (  I  ,  JA)**2. ) 

„ _ GO  TO  57  _ _ _ _ _ 

56  DPRE  aoi  i,  J)=0.0 

57  CONTINUE 

WR  I  T  E ( 6 , 1 30 ) _ 

130  FORMAT ( 1 -  * ,2 OX, * DIMENSIONLESS  PRESSURE  PROF  I •  , 

l'LE*  ) 

_ KRIT£I6«132I _ _ _ _ _ __ _  ... _ 

132  FOKMAT (»©• »2X, 1 TIME  * , 20X,  4  P ( X, T ) **2-PF**2 /P£**2-PF*» ’  , 

1*2) 

DO  58  1  =  1, IA _ _ _ 

5  8  WRITE(6f131)TU),  ( DPRES  S  (  I  ,  J  )  ,J=1,JA) 

131  FORMAT  (  '  •  ,F7.3,2X, 13(1X,F6.4)  ) 

_ . _ . PR=py..atJA) . ; . . _ . . . . . . . . — . . . . 

CALL  VISCN2tPR,TE,VN2i 
DO  6  5  1  =  1  ,  IA 

il  T  (  I  )  =  (  ,633E-03»P0(  1,  JA)  *PEkML*T  (!))/(  PH  I  »  VN  2  *X  (  J  A  )  **2 

1) 

DO  65  J= 1 , J  A 

_ DX(  J  )  =  XtJi/XUA) _ _ _ _ _ 

DP ( I  ,  J)=PL(  I  ,  J  )  /PL (  1 , JA) 

65  CONTINUE 

WRITEI6, 140) _ _ _ 

140  FORMAT ( 1 0 1  , 30X , ' P I  X , T ) /PF  »  ) 

WRITE (6 , 141 J ( OX ( J) , J= 1, JA) 

141  FORMAT (  IPX,  13<  1X.F6.4I) _ __ _ _  _ _ 

wR  I  T  E  (  6 , 1  <t3  ) 

143  FORMAT (IX, ' DT I  ME  *  ) 

DO  70  1  =  1  ,  I  A _ ; _ 

7  0  WRITE  (6, 142)  OH  I  )  ,  (DPI  1  ,  J)  ,  J=l,  JA) 

142  FORMAT!*  • ,F7.4,2X, 13 ( 1X,F6.4) ) 

60  CONTINUE . . . . . . __ . .  . 

STOP 

END 
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SUBROUTINE  PRCQEF ( PH I , T E , W ) 

DIMENSION  W ( 20  » 20 ) »  WPRS ( 3 ) 

COMMON  B(  lOl.CilO)  »0(  10)  15)»P  H  ( 13)  ,  T i  IS)  ,X(  15)  »  TP(  1 

10) » PERM (15), BrtB 
CGMMON  IG» JJ,JA,LL,LK 

JAM1=JA-1 _ 

DO  12  L=2,JAM1 
DO  12  J- 1 , J A 
. I  F { J . LT . L— 1 ) GO  TO  11 . 

IF( J.€T*i+llGiJ  fO  II 
IF(L,EU.J)GO  TO  7 

IFt J.GT.L iGU  TO  9 _ 

L  A  =  L- 1 
DO  2  1=1,3 

.  PR-P  <  LA  )  . . . . . . . . . . ___ . 

CALL  C0MPN2 I PR,TE, ZN2 ) 

CALL  VI  SC  N  2  ( P  R  ,  T  E  »  V  N2 ) 

WPRSI  I )  =l.i ( VN2*ZN2) _ 

L A=  LA+ 1 
2  CONTINUE 

_ Li=a _ _ _ ; . . . . . - . - . - _ - . . . 

CALL  PERMEF ( L , I  I , EFFK ) 

DELwPS=WPRS( 3)— WPRS ( 1 ) 

PJP11— EFFK » DEL WPS _ 

P JP 13=4 . *WPRS  C  2 ) *EFFK 
I  1  =  2 

CALL  PERMEFiL,  1 1  »EFFK) _ _ _ _ _ _ _ 

PJP12=WPRS(2)*EFFK*2. 

I  F { J J  *  EQ • 2  ) GO  TO  fa 

W(L.J)=PJP11+PJP12-PJP13 _ 

GO  TO  12 

6  W(L , J)=-PJP11-PJP12+PJP13 

GO  TO  12 . .  . 

7  D£LXSQ=X4L+1  >*X(L)-X(L  )  * X  (  L  I -X  (  L  +  1 )  * X  (  L - 1 )  +  X  (  LI*X(L-1) 
deLt=t (LK+1 )-T l LK) 

PK=P ( L  ) _ 

CALL  C0MPN2I PR,TE,ZN2) 

PJ  1=.  12fa4E+05*PHI*DELXSQ/  (P(  JA)  *DELT*ZN2) 

P J2=2 . *PJP1 3  . . . 


I  F { J J  »  EO • 2 ) GO  TO  8 
W(L  ,  J)=PJ1  +  PJ2 

GO  TO  12 


8 

9 

W(L, J)=PJ1-PJ2 

GO  TO  12 

1F{  JJ.Ey,.2)G0  TO  10 

10 

W  (  L  ,  J  )  =-(  PJPU+PJP12  +  PJP13) 

GO  TO  12 

W  (  L  ,  J  )  ■=  P  JP1  1+PJP12+P  JP13 

GO  TO  12 

11 

W(L,J) —0»0 

12 

CONTINUE 

RETURN 

END 
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S16RGUT I NE  CUMPN2 ( PX , T£ , ZN2 ) 

DIMENSION  ALIPC 10) , U  10 ) 

_  COMMON  BC10i«C(10) «D( 10) »P< 15KPMC15) «TC 15) «X( 131 » TP ( 1 
10) f FERMI  lb) ,Bab 
COMMON  I G , J J  ,  J A  » LL , LK 

P A=PX7 14. 696 _ 

ZN2  =  0 . 0 

T EK- { TE-32 . ) *(5./9.  )+273. 1 

. . . DO . lb . 1=1, 1G . . . . . . 

ALIPC I  )  =  1.0 
DO  10  J  =  1  ,  IG 

IF(I.£Q.J)GG  TO  10 _ 

ALIP  (I)  =  ALIPm  *(  (  T  EK-TP  (  J  )  )  /  (  TP(  I  )-TPC  J  )  )  ) 

10  CONTINUE 

1  {  1  )  =  1  *±fi  (  I  1  *PA+C  (  I  ) *PA**2.+D(  1  )  *PA**3-  _ ______ _ 

IFITEK.EQ.TPC U )  GO  TO  16 

15  ZN2=ZN2+ALIP< I )*Z(I) 

GO  TO  17 _ 

16  ZN2=Z(I) 

17  RETURN 

_JEljQ . . . ; _ _ _ _ _ _ _ — - - - - - - - 
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SUBROUTINE  V  I SGN2 I PX , TE , VN2 ) 

COMMON  B  (  IQ) » C I  10)  »D(  10) » P  {  15) »PM( 15)  »T (15) »  X  (  1 5  ) »  T  P  (  1 

_ _ _ _ _ _ _ _ _ _ 

COMMON  1G» JJ » JA,LL,LK 
P  A  =  P  X/ 1 4 . 696 

PMC  =  P  A—  1  » _ 

UR  AT  10^=1  .+.895E-03*PMC+.612E-06*PMC**2.+3.997E-08*PMC* 
1*3. 

VI  SN2-=Uft AT  1 0*1 «  778E-02  _ _ _ _ _ _ _ _ 

OELT  EC= ( T  E-32 . ) *(5./9.)-25. 

VN2=VISN2+0ELTEC*.455E-C4 

RETURN 


END 


<  fit  >  .■  S  ,  v 
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SUBROUTINE  P EKMEE { L , I  1 , EFEK ) 

COMMON  81  10)  ,C{  10)  ,□(  10)  ,P(  15  i  ,PM(  15)  ,T(  15)  ,X(  15)  ,TP{  1 

IQ)  ,P ERMt 15} «  88B  . . . . . . . 

COMMON  1G» JJf JA,LL,LK 
IF ( I  I . EQ . 2 ) GO  TO  5 

OELXI=( X( L )  -  X  {  L  - 1 ))/2. _ 

DELX2=<XiL*li-X(L) )/2, 

EFFK— ( QELXl^PERMI L- 1 ) +DELX2*PERM( L) ) / ( 0ELX1 +0ELX2 ) 

GO  TO  10 _ _ _ _ _ _ _ 

5  EFFK=PERM (L)-PERM(L-l ) 

10  RETURN 
END 
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•  -M  -  •  -t  >•  • 
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SUBRuUT I Nt  SQRGCT {  A  ,  B  A  ,  Y  ) 

DIMENSION  A(20#20) »AA( 20,20 ) fAL( 20,20 ), AM (2 0,20) ,AU( 20 
1«20I,AV(20,  gQJ.  ,  BAi  151,881  15  J  ,BC (20  »2QJ..J3JLE115-i  tPCENT  (  i 
15) ,  Y( 15) 

COMMON  6( 10) ,C(10) ,D{ 10) ,P( 15) ,PM( 15) ,T( 15) ,X( 15) ,TP( i 
_ 10)  ,PERM( 15) ,  BBB _ 

COMMON  iG , JJ » JAtLL ,LK 

151  FORMAT (1H  ,29H  THE  METHOD  IS  NOT  APPLICABLE) 

— _ M=l.L. . . . . . . . . . . . . . . - . - . 

DC  1  1  =  1 »  M 
DO  2  J^L  ,  M 

2  AA(I»J)-=A(1,J) _ 

1  68(1 )  =  BA (  I  ) 

C 

C  CALCULATION  OF  L  AND . U _ _ _ _ _ _ _ 

DO  10  5  K  =  2  t  M 
KK=K-1 

_ DO  104  J- 1  i.KH _ 

AL ( J  *  K ) —  0  •  0 
AMI J,K)=0.0 

. A.U,(.K..t..^l..=  Q.*.IL_ _ _ _ _ _ _ _ _ _ _  _ _ — . 

1 C 4  AV(K,J)=0.0 
1C5  CONTINUE 

_ IF  t A( 1, 1 J .Ew.D.  )  GO  TO  1.5.Q _ 

106  X  I  S—  A  B  S  (  A  ( 1,1)  ) 

AL ( 1 , 1 )=SQRT (XIS) 

_ AU  (  1,1)  =  A  ( 1, 1 )  / AL  (1 , 1) _ _ _ _ _ . _ _ - 

DO  107  J=2,M 

AU ( 1 , J ) = A ( 1, J)/AL(1,1 ) 

1 C 7  AL(Jtl)  — A(Jt 1)/AU(1,1) _ 

DO  115  N— 2  7 M 
KN=K~1 

VALU£  =  Q  .  Q . . . . . . . . . . . . . _ _ . . . 

DO  108  J— 1 , K  K 

1  C  8  VALLE=VALUE-»-AL(K,  J)*AU(  J,K) 

_ ZX=A ( K  «  Ki -VALUE _ 

IEIZX.EQ.O.)  GO  TO  150 
1C9  ZXX— AB  S I ZX ) 

_ _ _  AL(K,K)=SQRT(ZXX) _ _ _ _ _ . _ _ _ _ _ _ _ _ 

AU ( K , K ) =  Z  X/ AL ( K , K ) 

KP-=  K+  1 

_ KK^=K—  1  _ _ 

IFCKP.GT.M  )  GO  TO  115 
DO  112  i=KP,M 

zv=o.o . . 

ZW=0.0 

DO  110  L P  = 1 , K K 

_ ZV=ZV+AL ( K,LP )*AU( LP,  I  ) _ 

110  ZW=ZW+Al( IfLP)*AUiLP,K) 

AU(K, 1 )=i A(K, I )-ZV)/AL(K,K) 

1 12  AL i 1 >K)  —  (AC  I )-ZW) / AU ( K , K ) . . 

115  CONTINUE 
C 

C_ L  AND  U  ARE  CALCULATED 
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C 

PROCEEDING  TO  CALCULATE  L-INVERSE 
DC  119  K=  1 ,  M 

AM(K,K)=1 .0/AL{K,K) 

AND  U- INVERSE 

1 19 

AVI K,K)=1.0/AU(K,K> 

DO  125  K=2,M 

KK=K~1 

- 

DO  122  J=  1 , KK 

ZQ=0.0 

DO  120  L= J  » KK 

120 

122 

125 

ZQ=ZQ  +  AU  K,L)*AN(L, J) 

AM { K i J ) ZQ/ AL ( K t  K) 

CONTINUE 

I JN=M  + 1 

DO  135  KL  =  1 , IJN 

K=I JN-KL  . . . . . 

126 

IF(K.LE.l)  GO  TG  136 

DO  130  J  K=  1 ,  K 

J— K— JK 

127 

I  F ( J  .  L  T  .  1 )  GO  TO  135 

ZR=Q . 0 

JP=J+1 

128 

1  30 

DO  128  l=  JP , K 

ZK=ZR+AU(JfL)*AV« L,K) 

AVI J,K)=-ZR/AU(Ji J) 

C 

C 

135 

CONTINUE 

4 

PROCEEDING  TO  CALCULATE  G-INVERSE 

=  U- INVERSES- INVERSE 

736 

DO  140  K=1*H 

DO  140  J=  1 »  M 

A ( J , K ) =0*0 

137 

140 

DO  137  L=1,M 

A I  J  »  KJ  =  A I J,KI+AVi J,L) *  A  M I  L »  K  ) 
CONTINUE 

150 

152 

GO  TC  152 
kRITE ( 6, 151 ) 

CONTINUE 

DO  153  I-1»M 

Ylli=0. 

DO  153  J=  1 »  N  .  _  .  „ _ _ _ _ _ _ _  _ - 

153 

Y(  I  )  =  Y( I )+A( I , J3*8AI J ) 

DO  9  J=  1  ,  P 

BGB=0 . 

12 

DO  12  1  =  1,  M 

B0B=B08+AA{ J,IJ*Y(I) 

CIFI J)=bB( J)-6G8 

9 

202 

PCENTI J)=QIFI  J) *100. /BBC J ) 
CONTINUE 

RETURN 

END 
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APPENDIX  D 


GRAPHICAL  RESULTS 


p(x,  t),  (PSIA) 
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FIGURE  D-l 

P  R  E  SSUR  E  DI STRIB U  TI  ON 
BEREA  1:  CONSTANT  TERMINAL  PRESSURE  WITH 
CONSTANT  PRESSURE  AT  EXTERNAL  BOUNDARY 


x/L.  (DIMENSIONLESS) 
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FIGURE  D-2 

PRESSURE  DISTRIBUTION 
I’) ERE  \  1 :  CONSTANT  TERMINAL  PRESSURE 
WITH  SEALED  EXTERNAL  BOUNDARY 


x/L.  (DIMENSIONLESS) 


p(x,t),  (PSIA) 
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FIGURE  D-3 

PR  ESSURE  DioT  R  IB  UT I  ON 

BEREA  1A:  CONSTANT  TERMINAL  PRESSURE  WITH 
CONSTANT  PRESSURE  AT  E  XTERNAL  BOUND  ARY 


x/L,  (DIMENSIONLESS) 
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FIGURE  D-4 

PRESSURE  DISTRIBUTION 

BEREA  1A:  CONSTANT  TERMINAL  PRESSURE  WITH 
SEALED  EXTERNAL  BOUNDARY 


x/L,  (DIMENSIONLESS) 


P(x,t),  (PSIA) 
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FIGURE  D-5 

PRESSURE  DISTRIBUTION 
BEREA  2  :  CONSTANT  TERMINAL  PRESSURE  WITH 
CONSTANT  PRESSURE  AT  EXTERNAL  BOUNDARY 


x/L,  (DIMENSIONLESS) 


(VIScl)  ‘(}‘x)d 
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FIGURE  D  6 

PRE  SSURE  DISTRIBUTION 
BEREA  2:  CONSTANT  TERMINAL  PRESSURE  WITH 
SEALED  EXTERNAL  BOUNDARY 


x/ 1 (DIMENSIONLESS) 


P(x,t),  (PSIA) 
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FIGURE  D-7 

PRESSURE  DISTRIBUTION 
SERIES  CORE:  ALUNDUM-BERE A 
CONSTANT  TERMINAL  PRESSURE  WITH  SEALED 
EXTERNAL  BOUNDARY 


x/L,  (DIMENSIONLESS) 
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FIGURE  D-8 

PRESSURE  DISTRIBUTION 
LIMESTONE  CONSTANT  TERMINAL  PRESSURE  WITH 
CONSTANT  PRESSURE  AT  EXTERNAL  BOUNDARY 


x/L,  (DIMENSIONLESS) 


(VISd)  *ulx)d 
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FIGURE  D-9 

PRESSURE  DISTRIBUTION 
BEREA  1:  CONSTANT  TERMINAL  PRESSURE  WITH 
CONSTANT  PRESSURE  AT  EXTERNAL  BOUNDARY 


0.2  0.4  0.6  0.8  1.0 


x/L,  (DIMENSIONLESS) 


p(x,t),  (PSIA) 
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FIGURE  D--10 
PRESSURE  DISTRIBUTION 
SERIES  CORE:  ALUNDUM  BEREA 
CONSTANT  TERMINAL  PRESSURE  WITH  CONSTANT 
PRESSURE  AT  EXTERNAL  BOUNDARY 


x/L,  (DIMENSIONLESS) 
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